1. IntroductionIn medical applications, computed tomography (CT) has been widely used to reconstruct the detailed structure of objects. A CT imaging system can be modeled as:
where
A denotes the system weight matrix,
p denotes the projection data collected by the detector, and
v denotes the linear attenuation coefficients of the reconstructed object. The task of CT imaging is to reconstruct
v from
p, x-ray computed tomography (CT) uses attenuation coefficients of an object to reconstruct interior structures. However, the attenuation values change little between different biological soft tissues, so tissue characteristics are indiscernible in x-ray CT images. In-line phase-contrast (IL-PC) imaging is a new imaging method that generates high spatial resolution and super contrast of weakly absorbing objects compared with conventional attenuation-contrast imaging.
[1] In practical biomedical applications, IL-PC is combined with CT, and IL-PC-CT can reconstruct images of soft samples. IL-PC-CT images, reconstructed from the two-dimensional (2D) x-ray projection data acquired with a flat-panel detector, are often corrupted by serious ring artifacts, because of malfunctioning and miscalibration of elements in the x-ray detectors as well as impurities in the scintillator crystals.
[2] As visualization and quantification of the anatomic and pathological features are influenced by these ring artifacts, removal or reduction of these artifacts is essential. Ring artifact removal or reduction without impairing the image quality is still a challenging problem for the researchers. To overcome the challenge, several methods have been reported so far to correct the ring artifacts from CT images. The flat field correction is a common method to reduce ring artifacts.
[3] However, ring artifacts are not removed completely by this method when different camera elements have intensity dependency, nonlinear response functions or the incident beam has time-dependent non-uniformities.
[2] Another method for ring artifacts suppression is based on moving the sample or the detector system during the acquisition in defined horizontal and vertical steps.
[4,5] However, this requires very precise translation motor components. The third method for correcting ring artifacts is the image-processing based on sinogram preprocessing and reconstructed images postprocessing. Some reported methods dealing with the raw sinogram aim at detecting and correcting the spurious lines in the sinogram before applying the reconstruction process.
[6,7] The post-processing approach may be tempting to detect and remove the ring artifact in the image domain rather than the sinogram domain.
[8,9]TV regularization has been a very active research topic in CT reconstruction.[10–13] One example in the field of image reconstruction is the isotropic TV norm.[14,15] Later, the isotropic TV minimization approach was addressed in in-line phase-contrast x-ray tomography.[16] Although the isotropic TV norm was proven effective for low-dose CT and interior tomography, the isotropic TV minimization is unfit for limited-angle CT. A weighted anisotropic TV minimization method[17] is introduced to address the limited-angle CT problem using a simultaneous algebraic reconstruction technique (SART) algorithm[18] and the gradient descent method.[19]
The TV model falls into two categories: isotropic TV and anisotropic TV.[20,21] The mathematical definitions for both isotropic and anisotropic TV are given in the discrete setting. v is denoted as the column vector by a lexicographical ordering of a two-dimensional (2D) image. In our study, the isotropic and anisotropic TV are respectively defined as[22]
where
is the first-order finite difference operator in the horizontal and vertical directions.
The TV term (isotropic or anisotropic) is essentially the
-norm of the gradient of the image.[23] Thus, the sparsity constraint
-norm forces the ring artifacts variables to have only a few not null components, since the
-norm is a good approximation of the sparsity-inducing
-norm. Similarly, in this study, the ring artifacts can be wide and faint, which can be corrected by the TV term. The isotropic TV equals to the square root of the first-order finite difference operator in the horizontal and vertical directions, which is similar to the
-norm of a complex vector.[24] Thus, first-order finite difference operator in the horizontal and vertical directions are unseparated, which leads to the mutual interference. When the images become more complex, i.e., images that contain complex textures and gradual intensity like liver images, isotropic total variation can often produce cartoon-like approximations; for the result see Subsection 3.1. The main idea of anisotropic TV is to separate the first-order finite difference operator in the horizontal and vertical directions to avoid the interference between them, hence the anisotropic total variation regularization can preserve images sharp with corners,[25,26] and natural images usually contain a lot of edges and other anisotropic features, since the anisotropic TV is a good sparsity approximation of natural images. In order to reduce visible ring artifacts while preserve the image section details(edges, corners, and other anisotropic features), based on the advantages of anisotropic TV outlined above, we propose to correct the ring artifacts by use of an anisotropic TV term during the image reconstruction process; for the result see Subsection 3.2.
There are two differences between our method and the weighted anisotropic TV minimization method.[17] Firstly, the weighted anisotropic TV minimization method is introduced to address the limited-angle CT problem, while our aim is to correct image ring artifacts and reconstruct an image from the complete-view without introducing weighting factors for anisotropic TV minimization. Secondly, unlike the calculation method of the weighted anisotropic TV minimization,[17] in our study, the anisotropic TV minimization method is introduced to correct image ring artifacts using a simultaneous algebraic reconstruction technique (SART) algorithm[18] and the alternating direction method of multipliers (ADMM).[27]
The main features and contributions of this work are summarized as follows.
The rest of this paper is organized as follows. In Section 2, we will describe an anisotropic TV-L2 model and its implementation based on the ADMM.[27] In Section 3, the real liver section synchrotron data results are presented for the anisotropic TV-L2 model and demonstrate the effectiveness of the proposed method for correction of ring artifacts. Finally, concluding remarks are provided in Section 4.
2. Method2.1. In-line phase-contrast imagingThe experiments were performed in BL13W1 of the Shanghai synchrotron radiation facility (SSRF). A liver sample was placed on the shelf, and its image was reconstructed using an IL-PC x-ray imaging setup. Figure 1 shows the schematic of the IL-PC x-ray imaging setup.[28] The monochromator crystal was a combination of an Si (111) orientation crystal and an Si (311) orientation crystal. A high-precision sample platform was used to position the sample and rotate the sample axis perpendicular to the beam. The x-ray beam energy in the experiments was set at 20 keV, and the detector employed an x-ray CCD camera system with
per pixel. The projections were recorded with a sample-to-detector distance of 1 m and exposure time of 10 ms (as shown in Fig. 2(a)). In addition, 20 flat field images and 10 dark field images were collected.
Figure 2(a) shows a sinogram image of the projection data from the liver section. The data were collected by taking 1024 projections with total 180° angle of view. In Fig. 2(b), the magnified ROI of the sinogram reveals vertical lines that are responsible for ring artifacts. In this study, our aim is to correct ring artifacts during the reconstruction process.
2.2. Algebraic reconstruction methodTo solve the problem in Eq. (1), we use SART to enforce consistency with the projection data. The formula of SART can be expressed as[18]
where
K indicates the iteration number,
means the
j-th 2D pixel,
is the
n-th 2D pixel contribution along the
i-th ray,
is the measured projection, and
is the relaxation parameter.
2.3. Anisotropic TV-L2 modelIn our study, an image is reconstructed as a solution for the following anisotropic TV-L2 model:
where
is the regularization parameter.
A new reconstruction algorithm is developed based on proximal forward–backward splitting (PFBS),[29] which completely decouples data fidelity minimization and image regularization problems. As a result, a CT image reconstruction algorithm step without any regularization operators can be derived for data fidelity minimization, and then sole image regularization or denoising subproblems without projection operators can be rapidly solved using ADMM.[27] The detail of the PFBS algorithm regarding how to decouple data fidelity minimization and image regularization problems was introduced by Guo in 2016.[30] According to the PFBS optimization algorithm, the anisotropic TV-L2 image reconstruction model mainly consists of the following two iterative steps with the initial image
:
In Eq. (
6),
v is the image to be reconstructed,
p is the projection data,
A the system matrix or a projection operator discretized from the x-ray transform, the step for computing
is the gradient descent update with a step size of
for the problem
. There are many methods to solve the linear inverse problem for
v iteratively. In this study, the SART algorithm is used to solve the subproblem Eq. (
6).
For the sake of convenience, we can rewrite the subproblem Eq. (7) as
where
and
.
Next, we provide a solution algorithm for the subproblem Eq. (8) with anisotropic TV regularization Eq. (3) based on ADMM. We adopt a splitting method from Ref. [31] by introducing two auxiliary variables y = f, z = Df, and the optimization is rewritten as
Its corresponding augmented Lagrangian function is
where
λ is the regularization parameter, vectors
μ and
ρ are the Lagrange multipliers, and
δ and
θ are the penalty parameters. We use the alternating direction method to iteratively solve the following subproblems
The Lagrange multipliers are expressed as
The subproblem is to minimize the augmented Lagrangian function; these subproblems are investigated one by one.
In summary, the implementation steps of the TV-L2 model based on SART-anisotropic TV are given as follows:
2.4. Comparison method2.4.1. SART based on isotropic TV (SART-isotropic TV)The CT reconstruction problem given by Eq. (1) can be solved by the compressive sensing-based reconstruction method to minimize the image isotropic TV regularized by the projections. This process is equivalent to solving the following optimization program:
The model described by Eq. (
22) was investigated using the alternating minimization procedure from:
[14]The gradient descent is controlled via the step size parameter β. The SART step loop is executed K times.
denotes the total number of gradient descents.
2.4.2. SART based on sinogram preprocessing with the wavelet-FFT filter technique(SART-wavelet-FFT filter)In order to correct ring artifacts, a fast, powerful and stable filter based on combined wavelet and FFT transforms was presented to preprocess the sinogram.[7] The steps of the filter based on combined wavelet and FFT transforms are as follows:
Step 1 in the first step, a tight condensation of the stripe information by combining the wavelet and FFT transforms.
Step 2 it then performs a damping of the relevant coefficients by using a Gaussian function.
When the sinogram preprocessing is complete, the output sinogram can be used to reconstruct the image by the SART algorithm. The reconstruction algorithm based on sinogram preprocessing is referenced as the SART-wavelet-FFT filter.
3. Numerical resultsThe TV-L2 model based on SART-anisotropic TV can be used to correct the ring artifacts for IL-PC-CT. We compare the algorithm with closely related methods: SART, SART-isotropic TV, and SART-wavelet-FFT filter. A real sample was used to evaluate our algorithm, and all the algorithms were implemented in C++ language and matlab language on a PC (4.0-GB memory, 3.4-GHz CPU). Figure 1 shows a schematic diagram of an in-line phase contrast x-ray imaging setup. The projection data were collected from evenly distributed over pi-view in the BL13W1 of the Shanghai Synchrotron Radiation Facility (SSRF) in China. The images with a size of 1024 × 1024 pixels were reconstructed by SART, SART-isotropic TV, SART-wavelet-FFT filter, and SART-anisotropic TV. The relaxation parameter for SART was set to
. The SART algorithm performed iterations K = 10.
3.1. The ring artifact correction method based on SART-isotropic TVThe TV-L2 model based on SART-isotropic TV in Eq. (22). We will illustrate some properties of the SART-isotropic TV based on numerical studies. The step size parameter β and the total number of gradient descents
are chosen in order to obtain visually satisfactory results. We now discuss the impact of β for SART-isotropic TV. Particularly, we choose β from
when we fix the total number of gradient descents.
As for SART-isotropic TV, the impact of β is illustrated in Fig. 3, in which we fix the total number of gradient descents
. For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 3(d)–3(f), which shows that SART-isotropic TV is sensitive to the step size parameter β. We can see that large β leads to over-smoothed, and hence a moderate
should be chosen for the best visualization.
We now discuss the impact of
for SART-isotropic TV. Particularly, we choose
from
when we fix the step size parameter β = 0.02.
As for SART-isotropic TV, the impact of
is illustrated in Fig. 4, in which we fix the β = 0.02. For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 4(d)–4(f), which shows that the ring artifact was suppressed as the increase of iteration index, and hence the
should be chosen for the best visualization.
3.2. The ring artifact correction method based on SART-anisotropic TVThe TV-L2 model based on SART-anisotropic TV in Eq. (5). We will illustrate some properties of the SART-anisotropic TV based on numerical studies. SART-anisotropic TV is controlled via the regularization parameter λ. We choose different regularization parameter
,
,
when we fix the maximum iteration number k = 110 of the solution of subproblems (11)–(15).
As for SART-anisotropic TV, the impact of λ is illustrated in Fig. 5, in which we fix the k = 110. For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 5(d)–5(f). Figure 5(e) shows that
should be chosen for the balance between the ring artifacts reducing and the liver section details preserving.
We now discuss the impact of maximum iteration number k of the solution of subproblems (11)–(15) for SART-anisotropic TV. Particularly, we choose k from {110,140,170} when we fix the regularization parameter
. As for SART-anisotropic TV, the impact of k is illustrated in Fig. 6.
For the sake of better visualization, one region of interest (ROI) was magnified in Figs. 6(d)–6(f), which shows that the SART-anisotropic TV succeeds in removing almost all visible rings to an appreciable level as the increase of iteration index. Figure 6(f) shows that the image is visually more satisfying, and hence the k = 170 should be chosen for the best visualization.
3.3. Comparison of ring artifact correctionFirst, we compare SART-anisotropic TV with two methods: SART and SART-isotropic TV. Figures 7(a)–7(c) with a size of 1024 × 1024 pixels were reconstructed by different methods using 1024 projections. Figure 7(b) was reconstructed by SART-isotropic TV with β = 0.02 and
. Figure 7(c) was reconstructed by SART-anisotropic TV with
and k = 170.
Figures 7(d)–7(f) illustrate this effect by means of the difference images between images reconstructed by SART and methods: SART, SART-isotropic TV, and SART-anisotropic TV. Figure 7(e) shows that liver structures are visible in the difference image, which implies that SART-isotropic TV cannot preserve the liver section details. As can be seen from Fig. 7(f), no liver structures are visible in the difference image, which implies that SART-anisotropic TV can preserve spatial resolution.
Next, in order to provide further validation of the SART-anisotropic TV reconstruction accuracy, one region of interest (ROI) was magnified to allow a better detailed comparison. The magnified regions of the SART are illustrated in Fig. 8(d). The experiment result in Fig. 8(d) shows that the image reconstructed by SART is influenced by the ring artifacts. Figure 8(e) shows that ring artifacts correction conflicts with detail preservation for SART-isotropic TV. However, the results of reconstructions for the liver sample in Fig. 8(f) demonstrate that SART-anisotropic TV is effective and efficient for structure characteristics preservation on the ring artifacts correction process.
Finally, in order to make a comparison study among the three algorithms: SART, SART-wavelet-FFT filter, and SART-anisotropic TV, two regions of interest(I and II)were magnified, which are shown in Fig. 9(a).
As shown in Figs. 9(d) and 9(g), without sinogram preprocessing, the image reconstructed by SART is so strongly contaminated by ring artifacts, that in order to correct ring artifacts, the input sinogram is preprocessed using a wavelet-FFT filter. The wavelet-FFT filter algorithm,[7] in addition to the input sinogram, employs three parameters for the filtering process: the highest decomposition level L, the wavelet type, and the damping factor σ. The choice of these parameters permits the removal of most ring artifacts from the reconstructed slice, with minimal distortion at the image center. Higher order wavelets have a stronger effect on the center of the image, while decomposition levels smaller than 5 are instead not enough to remove all ring artifacts.[7] Taking this factor into account, the results have been obtained with the DB42 wavelet L =5 and σ = 2.4 (Fig. 9(b)). The experiment result in Fig. 9(h) shows that most ring artifacts are reduced from the reconstructed slice. However, Fig. 9(h) also shows that the image reconstructed by the SART-wavelet-FFT filter still contains a lot of noise. The experiment result in Fig. 9(e) shows that there are minimal distortions at the reconstructed image center(solid arrow). Compared with SART-wavelet-FFT filter, SART-anisotropic TV still gives a good indication of ring artifacts correction, which were shown in Figs. 9(f) and 9(i).